##########################################
# COVARIATE CHECK: EXCLUSION RESTRICTION #
##########################################
# Author: Kasia Nalewajko
# First created: 09 March 2023
# Replicated: 09 June 2024

rm(list = ls())

# LOAD PACKAGES -----------------------------------------------------------

if (!require("dplyr")) install.packages("dplyr")
if (!require("estimatr")) install.packages("estimatr")
if (!require("haven")) install.packages("haven")
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("stringr")) install.packages("stringr")

# LOAD DATA ---------------------------------------------------------------

load("./00 SUBMITTED/00 APSR final/04 replication_files/01 data/main.Rda")

# RUN THE MODELS -------

# Generate vector of outcomes
outcomes <- c("log(turnout_19_pct+1)",
              "log(empty_ballots_19_pct+1)",
              "log(nuance_extremegauche_19_pct+1)",
              "log(nuance_gauche_19_pct+1)",
              "log(nuance_centregauche_19_pct+1)",
              "log(nuance_centredroit_19_pct+1)",
              "log(nuance_droite_19_pct+1)",
              "log(nuance_extremedroite_19_pct+1)",
              "log(nuance_divers_19_pct+1)",
              "log(turnout_24_pct+1)",
              "log(nuance_extremegauche_24_pct+1)",
              "log(nuance_gauche_24_pct+1)",
              "log(nuance_centregauche_24_pct+1)",
              "log(nuance_centredroit_24_pct+1)",
              "log(nuance_droite_24_pct+1)",
              "log(nuance_divers_24_pct+1)",
              "log(turnout_32_pct+1)",
              "log(empty_ballots_32_pct+1)",
              "log(nuance_extremegauche_32_pct+1)",
              "log(nuance_gauche_32_pct+1)",
              "log(nuance_centregauche_32_pct+1)",
              "log(nuance_centredroit_32_pct+1)",
              "log(nuance_droite_32_pct+1)",
              "log(nuance_divers_32_pct+1)",
              "log(turnout_36_1_pct+1)",
              "log(empty_ballots_36_1_pct+1)",
              "log(nuance_extremegauche_36_1_pct+1)",
              "log(nuance_gauche_36_1_pct+1)",
              "log(nuance_centregauche_36_1_pct+1)",
              "log(nuance_centredroit_36_1_pct+1)",
              "log(nuance_droite_36_1_pct+1)",
              "log(nuance_extremedroite_36_1_pct+1)",
              "log(nuance_divers_36_1_pct+1)",
              "log(collabos_antijew1000+1)"
)

# Create empty containers
estimate <- rep(NA, 1*length(outcomes))
std_error <- rep(NA, 1*length(outcomes))
term <- rep(NA, 1*length(outcomes))
nobs <- rep(NA, 1*length(outcomes))
CI_lower <- rep(NA, 1*length(outcomes))
CI_upper <- rep(NA, 1*length(outcomes))
models_list <- list()
n <- length(outcomes)

# Loop over outcomes
for(i in 1:length(outcomes)){
  
  models_list[[i]] <- summary(lm_robust(as.formula(paste(outcomes[i],"~ log(mpf1000+1) + log(pop1911) + area_sqkm + longitude + longitudesq + latitude + latitudesq")),
                                        clusters = id_bureau,
                                        fixed_effects = as.factor(ar_name),
                                        data = main))
  
  estimate[[i]] <- models_list[[i]][["coefficients"]][1,1]
  std_error[[i]] <- models_list[[i]][["coefficients"]][1,2]
  term[[i]] <- models_list[[i]][["call"]][["formula"]][[2]][[3]]
  nobs[[i]] <- models_list[[i]][["nobs"]]
  CI_lower[[i]] <- models_list[[i]][["coefficients"]][1,5]
  CI_upper[[i]] <- models_list[[i]][["coefficients"]][1,6]
  
  print(i)
}


# PUT TOGETHER DATA FRAME FOR PLOTTING -------

models_df <- data.frame(
  outcome = outcomes,
  variables = as.character(term),
  estimate = estimate,
  std_error = std_error,
  nobs = nobs,
  CI_lower = CI_lower,
  CI_upper = CI_upper)

models_df$outcome <- stringr::str_remove(string = models_df$outcome, pattern = "^log\\(")
models_df$outcome <- stringr::str_remove(string = models_df$outcome, pattern = "\\+.*$")

models_df$outcome <- str_replace_all(string = models_df$outcome, pattern = "_19_pct", replacement = " (1919)")
models_df$outcome <- str_replace_all(string = models_df$outcome, pattern = "_24_pct", replacement = " (1924)")
models_df$outcome <- str_replace_all(string = models_df$outcome, pattern = "_32_pct", replacement = " (1932)")
models_df$outcome <- str_replace_all(string = models_df$outcome, pattern = "_36_1_pct", replacement = " (1936)")
models_df$outcome <- str_replace_all(string = models_df$outcome, pattern = "turnout", replacement = "Turnout")
models_df$outcome <- str_replace_all(string = models_df$outcome, pattern = "empty_ballots", replacement = "Casted empty ballots")
models_df$outcome <- str_replace_all(string = models_df$outcome, pattern = "nuance_extremegauche", replacement = "Extreme left votes")
models_df$outcome <- str_replace_all(string = models_df$outcome, pattern = "nuance_gauche", replacement = "Left votes")
models_df$outcome <- str_replace_all(string = models_df$outcome, pattern = "nuance_centregauche", replacement = "Centre left votes")
models_df$outcome <- str_replace_all(string = models_df$outcome, pattern = "nuance_centredroit", replacement = "Centre right votes")
models_df$outcome <- str_replace_all(string = models_df$outcome, pattern = "nuance_droite", replacement = "Right votes")
models_df$outcome <- str_replace_all(string = models_df$outcome, pattern = "nuance_extremedroite", replacement = "Extreme right votes")
models_df$outcome <- str_replace_all(string = models_df$outcome, pattern = "nuance_divers", replacement = "Other votes")
models_df$outcome <- str_replace_all(string = models_df$outcome, pattern = "collabos_antijew1000", replacement = "Anti-Jewish collaborators")

# PLOT ---------------------------------------------------------------

models_df %>%
  ggplot(aes(reorder(outcome, -estimate), 
             estimate
  )) +
  geom_hline(aes(yintercept = 0), linetype = 2, color = "gray20") +
  geom_linerange(aes(ymin = CI_lower, ymax = CI_upper),
                 linewidth = 1) +
  geom_point(size = 1.75
  ) +
  theme_bw() +
  theme(legend.title = element_blank(),
        panel.background = element_rect(fill = NA, color = "black")
  ) +
  labs(
    x = "Outcomes",
    y = "Estimates on WWI soldier death rates"
  ) +
  coord_flip()

# SAVE ---------------------------------------------------------------

ggsave(
  "A4_covtest_postWW1.png",
  plot = last_plot(),
  path = "./00 SUBMITTED/00 APSR final/02 supplementary_materials/figures",
  width = 20,
  height = 24,
  units = "cm",
  dpi = 300
)



